
% unit cost of investment
if (Cobb_Douglas)
        
    % the share parameters already sum to one
    share_g=a_g;
    share_m=a_m;
    share_f=a_f;
    share_y=a_y;
    
    pbar=(p./share_g).^share_g.*(P_c./share_y).^share_y.*(W_m./share_m).^share_m.* ...
         (W_f./share_f).^share_f;
     
    pbar_single=(p./share_g).^share_g.*(P_c./share_y).^share_y.*(W_m./share_m).^share_m;
    pbar(married==0)=pbar_single(married==0);

else

    % tau_m/g
    Phi_m=(a_g./a_m.*W_m./p).^(1/(rho-1));
    
    % tau_f/g
    Phi_f=(a_g./a_f.*W_f./p).^(1/(rho-1));
    
    % define something for single and married
    aux=a_m.*Phi_m.^rho+a_f.*Phi_f.^rho+a_g;
    aux_single=a_m.*Phi_m.^rho+a_g;
    aux(married==0)=aux_single(married==0);
    
    % Y_c/g
    Phi_c=(a_g.*a_h./a_y.*P_c./p.*aux.^(gma/rho-1)).^(1/(gma-1));
    
    % g/X
    Theta_g=(a_h.*aux.^(gma/rho)+a_y.*Phi_c.^gma).^(-1/gma);
    
    pbar=Theta_g.*(p+P_c.*Phi_c+W_m.*Phi_m+W_f.*Phi_f);
    pbar_single=Theta_g.*(p+P_c.*Phi_c+W_m.*Phi_m);
    pbar(married==0)=pbar_single(married==0);
    
    % expenditure shares
    share_m=Phi_m.*Theta_g.*W_m./pbar;
    share_f=Phi_f.*Theta_g.*W_f./pbar;
    share_g=Theta_g.*p./pbar;

    share_f(married==0)=0;
    share_y=1-share_m-share_f-share_g;

end

K=alpha.*beta.^(T-age+1).*delta2.^(T-age).*delta1;

% annual consumption, investment and leisure
if (price_reduction & reduce_W & t>=5)
    sum_cshare=(1-R.^(-(T_D-4)))/(1-1/R)+(1-delta2^(-(T-4)))/(1-1/delta2).*alpha.*(beta^(T-4)*delta2^(T-5)*delta1)+(1-R.^(-(T_R-4)))/(1-1/R).*(psi_m+psi_f);
else
    sum_cshare=(1-R.^(-T_D))/(1-1/R)+(1-delta2^(-T))/(1-1/delta2).*alpha.*(beta^T*delta2^(T-1)*delta1)+(1-R.^(-T_R))/(1-1/R).*(psi_m+psi_f);
end

c=1./sum_cshare.*W./p;
X=K./sum_cshare.*W./pbar;
l_m=psi_m./sum_cshare.*W./W_m;
l_f=psi_f./sum_cshare.*W./W_f;

% per child weekly investment
X=X/52;

% investment expenditure (weekly per child)
pbar_X=pbar.*X;
p_g=share_g.*pbar_X;
W_m_tau_m=share_m.*pbar_X;
W_f_tau_f=share_f.*pbar_X;
P_c_Y_c=share_y.*pbar_X;

% investmenet quantity
g=p_g./p;
tau_m=W_m_tau_m./W_m;
tau_f=W_f_tau_f./W_f;
Y_c=P_c_Y_c./P_c;

% weekly hours worked
hrs_m=100-tau_m-l_m/52;
hrs_f=100-tau_f-l_f/52;

% annual quasi difference in log achievement
qd_logPsi=delta1.*log(X.*52)+X_theta*phi_theta;

%% sample selection

if ((decomposition & ~baseline) | (price_reduction & ~quo))
    insample=insample_baseline;
else
    insample=(insample_data==1 & hrs_m>0 & (married==0 | (married==1 & hrs_f>0)));     

    % also check lifetime labor supply
    for i=1:N
        if (insample_data(i)==1)
            for t=T+1:T_R(i)
                x=expr_m0(i)+t;
                hrs=100-(psi_m(i)/sum_cshare(i)*W(i)/(W_m_base(i)*exp([x x^2]*param_W_m)))/52;
                if (hrs<0)
                    insample(i)=0;
                end
                
                if (married(i)==1)
                    x=expr_f0(i)+t;
                    hrs=100-(psi_f(i)/sum_cshare(i)*W(i)/(W_f_base(i)*exp([x x^2]*param_W_f)))/52;
                    if (hrs<0)
                        insample(i)=0;
                    end
                end
            end
        end
    end
    
end

qd_logPsi(insample==0)=nan;

p_g(insample==0)=nan;
W_m_tau_m(insample==0)=nan;
W_f_tau_f(insample==0)=nan;
P_c_Y_c(insample==0)=nan;

tau_m(insample==0)=nan;
tau_f(insample==0)=nan;
g(insample==0)=nan;
Y_c(insample==0)=nan;
X(insample==0)=nan;
pbar(insample==0)=nan;

%% model shares
share_m(insample==0)=nan;
share_f(insample==0)=nan;
share_g(insample==0)=nan;

for m=1:2
    for e=1:2
        
        eval(['share_model(m,e,1)=mean(share_m(ind' int2str(m) int2str(e) '),''omitnan'');']);
        if (m==2)
            eval(['share_model(m,e,2)=mean(share_f(ind' int2str(m) int2str(e) '),''omitnan'');']);
        end
        eval(['share_model(m,e,3)=mean(share_g(ind' int2str(m) int2str(e) '),''omitnan'');']);
    end
end

%% model levels
pbar_X(insample==0)=nan;
hrs_m(insample==0)=nan;
hrs_f(insample==0)=nan;

for m=1:2
    for e=1:2

        eval(['level_model(m,e,1)=mean(tau_m(ind' int2str(m) int2str(e) '),''omitnan'');']);        
        eval(['level_model(m,e,2)=mean(hrs_m(ind' int2str(m) int2str(e) '),''omitnan'');']);
        if (m==2)
            eval(['level_model(m,e,3)=mean(hrs_f(ind' int2str(m) int2str(e) '),''omitnan'');']);
        end
        
    end
end

%% detailed outcomes
if (decomposition)
    
    outcome=zeros(14,2);
    doutcome=zeros(14,2);
    
    for e=1:2
        j=0;
        
        for m=1:2
            
            % expenditure
            j=j+1;
            eval(['outcome(j,e)=mean(pbar_X(ind' int2str(m) int2str(e) ...
                '),''omitnan'');']);

            if e>=2
                doutcome(j,e)=(outcome(j,e)/outcome(j,1)-1)*100;
            end

            % price
            j=j+1;
            eval(['outcome(j,e)=mean(pbar(ind' int2str(m) int2str(e) ...
                '),''omitnan'');']);

            if e>=2
                doutcome(j,e)=(outcome(j,e)/outcome(j,1)-1)*100;
            end
            
            % total quantity            
            j=j+1;
            eval(['outcome(j,e)=mean(X(ind' int2str(m) int2str(e) ...
                '),''omitnan'');']);
            
            if e>=2
                doutcome(j,e)=(outcome(j,e)/outcome(j,1)-1)*100;
            end

            % mother's time
            j=j+1;
            eval(['outcome(j,e)=mean(tau_m(ind' int2str(m) int2str(e) ...
                '),''omitnan'');']);
            
            if e>=2
                doutcome(j,e)=(outcome(j,e)/outcome(j,1)-1)*100;
            end

            % mother's time expenditure
            j=j+1;
            eval(['outcome(j,e)=mean(W_m_tau_m(ind' int2str(m) int2str(e) ...
                '),''omitnan'');']);
            
            if e>=2
                doutcome(j,e)=(outcome(j,e)/outcome(j,1)-1)*100;
            end

            % mother's time expenditure share
            j=j+1;
            eval(['outcome(j,e)=mean(share_m(ind' int2str(m) int2str(e) ...
                '),''omitnan'');']);
            
            if e>=2
                doutcome(j,e)=(outcome(j,e)-outcome(j,1))*100;
            end
                        
            % quasi difference in log achievement
            j=j+1;
            eval(['outcome(j,e)=mean(qd_logPsi(ind' int2str(m) int2str(e) ...
                  '),''omitnan'');']);
            
            if e>=2
                doutcome(j,e)=(outcome(j,e)-outcome(j,1))*100;
            end
        end
    end
    
    % variance decomposition
    vars=zeros(3,4);
    
    % put back data for expenditure
    pbar_X=pbar_X_data;
    pbar_X(insample==0)=nan;
    X=pbar_X./pbar;
    
    vars(1,1:3)=[var(log(pbar_X),'omitnan') var(log(pbar),'omitnan') var(log(X),'omitnan')];
    vars(2,1:3)=[var(log(pbar_X(married==0)),'omitnan') var(log(pbar(married==0)),'omitnan') var(log(X(married==0)),'omitnan')];
    vars(3,1:3)=[var(log(pbar_X(married==1)),'omitnan') var(log(pbar(married==1)),'omitnan') var(log(X(married==1)),'omitnan')];

    for m=1:3
        vars(m,4)=(vars(m,1)-vars(m,2)-vars(m,3))/(2*sqrt(vars(m,2)*vars(m,3)));
    end
    
end

